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Abstract 

The purpose of this paper is to analyze and compute the early exercise boundary for a 



class of nonlinear Black-Scholes equations with a nonlinear volatility which can be a function 
of the second derivative of the option price itself. A motivation for studying the nonlinear 
Black-Scholes equation with a nonlinear volatility arises from option pricing models taking 
into account e.g. nontrivial transaction costs, investor's preferences, feedback and illiquid 
markets effects and risk from a volatile (unprotected) portfolio. We present a new method how 
to transform the free boundary problem for the early exercise boundary position into a solution 
^ of a time depending nonlinear parabolic equation defined on a fixed domain. We furthermore 

propose an iterative numerical scheme that can be used to find an approximation of the free 
boundary. We present results of numerical approximation of the early exercise boundary for 
various types of nonlinear Black-Scholes equations and we discuss dependence of the free 
boundary on various model parameters. 

d 
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1 Introduction 

H 

In the past years, the Black-Scholes equation for pricing derivatives has attracted a lot of attention 
from both theoretical as well as practical point of view. Recall that a European Call (Put) option 
is the right but not obligation to purchase (sell) an underlying asset at the expiration price E at the 
expiration time T. In an idealized financial market the price of an option can be computed from the 
well-known Black-Scholes equation derived by Black and Scholes in and, independently by 
Merton (see also EH . Dewynne et al. Q, Hull |16J). Assuming that the underlying asset follows 
a geometric Brownian motion one can derive a governing partial differential equation for the price 
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of an option. We remind ourselves that the equation for option's price V(S, t) is the following 
parabolic PDE: 

d t V + {r-q)Sd s V +\a 2 S 2 d 2 s V -rV = (1) 

where a is the volatility of the underlying asset price process, r > is the interest rate of a zero- 
coupon bond, q > is the dividend yield rate. A solution V = V(S, t) represents the price of an 
option at time t G [0, T] if the price of an underlying asset is S > 0. In this paper we shall focus 
our attention to the case when the diffusion coefficient a 2 may depend on the time T — t to expiry, 
the asset price 5* and the second derivative d 2 s V of the option price (referred to as V), i.e. 

a = a{S 2 d 2 s V,S,T-t). (2) 

A motivation for studying the nonlinear Black-Scholes equation (OQ) with a volatility a having 
a general form © arises from option pricing models taking into account nontrivial transaction 
costs, market feedbacks and/or risk from a volatile (unprotected) portfolio. Recall that the linear 
Black-Scholes equation with a constant has been derived under several restrictive assumptions 
like e.g. frictionless, liquid, complete markets, etc. We also recall that the linear Black-Scholes 
equation provides a perfectly replicated hedging portfolio. In recent years, some of these assump- 
tions have been relaxed in order to model, for instance, the presence of transaction costs (see e.g. 
Leland ll22l . Hoggard et al. ifTTl . Avellaneda and Paras B2l), feedback and illiquid market ef- 
fects due to large traders choosing given stock-trading strategies (Frey and Patie lfT2l . Frey and 
Stremme lfT3l . During et a/.jH, Schonbucher and Wilmott II2910 . imperfect replication and in- 
vestor's preferences (Barles and Soner O), risk from unprotected portfolio (Kratka [|20l . Jandacka 
and Sevcovic 020). One of the first nonlinear models is the so-called Leland model (c.f. 11221 ) 
for pricing Call and Put options under the presence of transaction costs. It has been generalized 
for more complex options by Hoggard, Whaley and Wilmott in IfTTl . In this model the volatil- 
ity a is given by a 2 (S 2 d 2 s V, S, r) = a 2 (l + Lesgn(<9fy)) where a > is a constant historical 
volatility of the underlying asset price process and Le > is the so-called Leland constant given 
by Le = y/2/nC / (ay/At) where C > is a constant round trip transaction cost per unit dollar of 
transaction in the assets market and At > is the time-lag between portfolio adjustments. 

If transaction costs are taken into account perfect replication of the contingent claim is no 
longer possible and further restrictions are needed in the model. By assuming that investor's pref- 
erences are characterized by an exponential utility function Barles and Soner (c.f. flU) derived a 
nonlinear Black-Scholes equation with the volatility a given by 

a 2 (S 2 8 2 s V, S, r) =a 2 (l + V(a 2 e rT S 2 8 2 S V)) (3) 

where ^ is a solution to the ODE: ^'(s) = (^(x) + l)/(2y / x^(x) - x), *(0) = 0, and a > is a 
given constant representing risk aversion. Notice that ^(x) = O(x^) for x —> and ^(x) = 0{x) 
for x — » oo. 

Another popular model has been derived for the case when the asset dynamics takes into ac- 
count the presence of feedback and illiquid market effects. Frey and Stremme (c.f. lfT3l WHi ) 
introduced directly the asset price dynamics in the case when a large trader chooses a given stock- 
trading strategy (see also [|29l ). The diffusion coefficient a is again nonconstant and it can be 
expressed as: 

a 2 (S 2 d 2 s V,S,r) = a 2 (l- gX(S)Sd 2 s V)- 2 (4) 
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where a 2 , g > are constants and X(S) is a strictly convex function, X(S) > 1. 

The last example of the Black-Scholes equation with a nonlinearly depending volatility is the 
so-called Risk Adjusted Pricing Methodology model proposed by Kratka in [20J and revisited by 
Jandacka and Sevcovic in [fT8l . In order to maintain (imperfect) replication of a portfolio by the 
delta hedge one has to make frequent portfolio adjustments leading to a substantial increase in 
transaction costs. On the other hand, rare portfolio adjustments may lead to an increase of the 
risk arising from a volatile (unprotected) portfolio. In the RAPM model the aim is to optimize the 
time-lag At between consecutive portfolio adjustments. By choosing At > in such way that the 
sum of the rate of transaction costs and the rate of a risk from unprotected portfolio is minimal one 
can find the optimal time lag At > (see lfT8l for details). In the RAPM model, it turns out that 
the volatility is again nonconstant and it has the following form: 

a 2 (S 2 d 2 V,S,r) = a 2 (l + ^(Sd 2 V)^ . (5) 

Here a 2 > is a constant and fi = 3(C 2 R/2tt)^ where C, R > are nonnegative constant repre- 
senting the transaction cost measure and the risk premium measure, resp. (see lfT8ll for details). 

Notice that all the above mentioned nonlinear models are consistent with the original Black- 
Scholes equation in the case the additional model parameters (e.g. Le, a, g, /S) are vanishing. If 
plain Call or Put vanilla options are concerned then the function V(S, t) is convex in 5 variable 
and therefore each of the above mentioned models has a diffusion coefficient strictly larger than a 2 
leading to a larger values of computed option prices. They can be therefore identified with higher 
Ask option prices, i.e. offers to sell an option. Furthermore, these models have been considered 
and analyzed mostly for European options, i.e. options can be exercised only at the maturity 
t = T. On the other hand, American options are more common in financial markets as they allow 
for exercising of an option anytime before the expiry T. In the case of an American Call option a 
solution to equation is defined on a time dependent domain < t < T, < S < Sf(t). It is 
subject to the boundary conditions 

V(0,t) = 0, V(S f (t),t) = S f (t)-E, d s V(S f (t),t) = l, (6) 

and terminal pay-off condition at expiry t = T 

V(S,T) =max(S-E,0) (7) 

where E > is a strike price (c.f. B71|2T1). One of important problems in this field is the analysis of 
the early exercise boundary Sf(t) and the optimal stopping time (an inverse function to Sf(t)) for 
American Call options on stocks paying a continuous dividend q > 0. However, an exact analytical 
expression for the free boundary profile is not even known for the case when the volatility a is 
constant. Many authors have investigated various approximation models leading to approximate 
expressions for valuing American Call and Put options: analytic approximations (Barone-Adesi 
and Whaley 0, Kuske and Keller JED, Dewynne et al. 0, Geske et al. flHE), MacMillan E51 - 
Mynemi ll26l ): methods of reduction to a nonlinear integral equation (Alobaidi [lj, Kwok j]2Tfl. 
Mallier et al. [|24l l25l . Sevcovic ll28l . Stamicar et al. Il3~0l ). Concerning numerical methods for 
solving free boundary problem we refer to the book by Kwok lETj and recent papers by Ehrhardt 
and Mickens [9] and Zhao et al. OTA . 
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The main goal of this paper is to propose a new iterative numerical algorithm for solving the 
free boundary problem for an American Call option in the case the volatility a may depend on 
the option and asset values as well as on the time T — t to expiry. The key idea of the proposed 
algorithm consists in transformation of the free boundary problem into a semilinear parabolic 
equation defined on a fixed spatial domain coupled with a nonlocal algebraic constraint equation 
for the free boundary position. Since the resulting parabolic equation contains a strong convective 
term we make use of the operator splitting method in order to overcome numerical difficulties. Full 
space-time discretization of the problem leads to a system of semi-linear algebraic equations that 
can be solved by an iterative procedure at each time level. 

The paper has the following plan: in the next section we transform the free boundary problem 
into a system consisting of a nonlinear parabolic equation defined on a fixed domain with time 
depending coefficients and an algebraic constraint equation for the free boundary position. In 
section 3 we present a numerical discretization scheme based on the idea of operator splitting. In 
the last section 4 we present several numerical results for nonlinear Black-Scholes equations with 
volatility functions a defined in ([3]) and [5]). We make a comparison to well-known methods in the 
case the volatility o is constant. Finally, we discuss dependence of the free boundary position with 
respect to various parameters entering expressions © and ©. 

2 Landau fixed domain transformation 

The main goal of this section is to perform a fixed domain transformation (referred to as Landau's 
transformation) of the free boundary problem for the nonlinear Black-Scholes equation (OQ) into a 
parabolic equation defined on a fixed spatial domain. For the sake of simplicity we will present a 
detailed derivation of an equation only for the case of an American Call option. Derivation of the 
corresponding equation for the American Put option is similar. 
Let us consider the following change of variables: 

r = T-t, x = In (g(r)/S) where g(r) = S f (T - r). 

Then r G (0,T) and x G (0, oo) iff S G (0, Sf(t)). The value x = corresponds to the free 
boundary position S = Sj(t) whereas x ~ +oo corresponds to the default value S = of the 
underlying asset. Following Stamicar et al. [l30l and Sevcovic ll28l we construct the so-called 
synthetic portfolio function II = H(x, r) defined as follows: 

U(x,t) = V(S,t)- Sd s V(S,t). (8) 

It corresponds to a synthetic portfolio consisting of one long positioned option and A = dgV 
underlying short stocks. Clearly, we have 

d x Tl = S 2 d 2 s V, d T Tl + ^d x U = -dt (V - Sd s V) 

Q 

where we have denoted g = dg/dr. Assuming sufficient smoothness of a solution V = V(S, t) to 
(OQ) we can deduce from (OQ) a parabolic equation for the synthetic portfolio function II = IL(x,t) 

d T U + (6(r) - \a 2 )d x Ii - \d x {a 2 d x U) + rll = 
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defined on a fixed domain x G R, t G (0, T), with a time-dependent coefficient 



and a diffusion coefficient given by: a 2 = a 2 (d x Il(x, r), ^(r)e~ x , r) depending on r, x and the 
gradient d x Tl of a solution II. Now the boundary conditions V(0, t) = 0, V(Sf(t),t) = Sj(t) — E 
and d s V{Sf{t), t) = limply 

IT(0,r) = -£, n(+oo,r) = 0, 0<r<T, (10) 

and, from the terminal pay-off diagram for V(S, T), we deduce 

n ( x, <,) = (-* *-* (11) 

I otherwise. 

In order to close up the system of equations that determines the value of a synthetic portfolio II 
we have to construct an equation for the free boundary position g(r). Indeed, both the coefficient 
b as well as the initial condition Il(x, 0) depend on the function g(r). Similarly as in the case of 
a constant volatility a (see 112811301 ) we proceed as follows: since Sf(t) — E = V(Sf(t),t) and 
d s V(S f (t),t) = 1 we have *S f (t) = d s V(S f (t),t)£S f (t) + d t V(S f (t),t) and sod t V{S,t) = 
along the free boundary S = Sf(t). Moreover, assuming 9 X II is continuous up to the boundary 
x = we obtain S 2 d 2 s V(S,t) -> d x U(0,r) and Sd s V(S,t) -> g(r) as S -> 5/(4)-. Now, 
by taking the limit S 1 — > S/(t)~ in the Black-Scholes equation dU) we obtain (r — q)g{r) + 
|a 2 <9n(0, r) - r(g(r) - £?) = 0. Therefore 

Q(t) = — + ^(^(0, r), e (r), r)^n(0, r) 
q zq 

for < r < T. Recall that in the linear case when a > is constant the initial position of the 
interface ^(0) is given by g(0) = rE/q if r > q > and g(0) = Eif0<r<q (see Dewynne et 
al. or Sevcovic [|28l ). We also recall that the value of ^(0) in the case r > q > can derived 
easily from the smoothness assumption made on d x Il at the origin x = 0, r = 0. Indeed, continuity 
of d x Tl at the origin (0, 0) implies lim T ^ + ^11(0, r) = 9^11(0, 0) = \im x _ t0 + d x Il(x, 0) = 
because IT(x, 0) = — E for x close to + . From the above equation for g(r) we deduce ^(0) = — 
by taking the limit r — > + . Henceforth, we restrict our attention to the case when the interest rate 
is greater then the dividend yield rate, i.e. 

0<g<r (12) 

leading to the initial position of the free boundary g(0) = rE/q. Putting all the above equations 
together we end up with a closed system of equations for II = H(x, r) and g = g(r) 

d T Tl + (6(r) - \a 2 )d x Yi - \d x {a 2 d x Yi) + rll = , 
n(0, r) = —E , II(+oo, r) = , x > , r G (0, T) , 

U(x0) = l~ E f ™ x <Hr/q) (13) 
v ' 1 \ otherwise, 
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where a = a(d x Il(x, r), g(r)e x , r) , b{r) = 44 + r — q and the free boundary position g 



q(t) 



T 



Sf(T — r) satisfies an implicit algebraic equation 



£ (r) = — + z d x IL{0, t) , with £ = — (14) 

q 2q q 

where r G (0, T). Notice that, in order to guarantee parabolicity of equation (TT3T) we have to 
assume that the function p i— > cr 2 (p, f?(T)e -a: , r)p is strictly increasing. More precisely, we shall 
assume that there exists a positive constant 7 > such that 

a 2 (p, £, r) + p<V 2 (p, £, t) > 7 > (15) 
for any £ > 0, r e (0, T) and p e R. 

Remark 1. In ll28l the author derived a single equation for the position of the free boundary g for 
the case when the volatility a = a is constant. In this case one can solve the initial-boundary value 
problem for the linear parabolic equation (fT3l) with spatially independent coefficients by means 
of one-sided sine and cosine Fourier transforms in the spatial x variable. The explicit formula 
for Il(x, r) together with equation (fl4l) enables us to conclude that the free boundary position g 
satisfies the following nonlinear weakly singular integral equation 

g(r) = r -^(l + ^= exp(-rr-(A T , + ln(r/g)) 2 /(2a 2 r)) 
q \ rV27rr 



+ 



a + — 1 — 



\ rE 1 /r — sJ -y/r — s 



ds (16) 



where the function A TS depends upon £> via A T)S = In g{f) — In g(s) + (r — q — \a 2 ) (r — s). 
The above integral equation can be solved by an iterative procedure based on a Banach fixed point 
argument (see ll28l for details). It is worthwhile noting that this approach cannot be applied to the 
case when the volatility may depend on the asset price S and/or the second derivative of the option 
price as the Fourier integral transform technique is no longer applicable. However, in section 
4 we shall use a solution computed from the nonlinear integral equation (fT6l) in order to make 
comparison of our iterative numerical scheme in the case when the volatility a is constant. 

Remark 2. We also present a formula for pricing American Call options based on the solution 

(II, g) to By © we have d s {S~ 1 V{S 1 1)) = -S~ 2 Tl (In (g(T - t)/S) ,T-t). Taking 

into account the boundary condition V(Sf(t),t) = Sf(t) — E and integrating the above equation 
from S to Sf(t) = g(T — t), we obtain the expression for the option price V(S, t): 



V(S, T-T) = -^Lr( g{r) -E+ [ * e x U(x, r)dx\ . (17) 

Qvn \ Jo J 

3 Discretization scheme. An iterative algorithm for approxi- 
mation of the early exercise boundary 

In this section we derive a full space-time discretization scheme for a numerical approximation 
of the problem (TT3T) . (TT4l) . Recall that in the case of a constant volatility there are, in principle, 
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two ways how to solve numerically the free boundary problem for the value of an American Call 
resp. Put option and the position of the early exercise boundary. The first class of algorithms is 
based on reformulation of the problem in terms of a variational inequality (see Kwok [|2T1l and 
references therein). The variational inequality can be then solved numerically by the so-called 
Projected Super Over Relaxation method (PSOR for short). An advantage of this method is that it 
gives us immediately the value of a solution. A disadvantage is that one has to solve large systems 
of linear equations iteratively taking into account the obstacle for a solution, and, secondly, the 
free boundary position should be deduced from the solution a posteriori. Moreover, the PSOR 
method is not directly applicable for solving the problem (Q])-© when the diffusion coefficient a 
may depend on the second derivative of a solution itself. The second class of methods is based on 
derivation of a nonlinear integral equation for the position of the free boundary without the need 
of knowing the option price itself (see e.g. Kuske-Keller lfT9l . Mallier and Alobaidi [[Tl l24l l25l . 
Sevcovic et al. Il28ll30l . In this approach an advantage is that only a single equation for the free 
boundary has to be solved; a disadvantage is that the method is based on integral transformation 
techniques and therefore the assumption a is constant is crucial. 

In our approach we make an attempt to take advantages of both above mentioned methods. As 
it was mentioned in the previous section we are going to solve the system of nonlinear equations 
(fT3l) with constraint (fl4l) . Because the volatility o may be nonconstant we cannot use integral 
transformation techniques in order to derive a single integral equation for g(r). However, the form 
of the system (TT3T) . (fT4j) allows for an efficient and fast numerical algorithm for computing of the 
early exercise boundary position g(r) = Sf(T — r). 

The idea of the iterative numerical algorithm for solving the problem (fT3l) . CHI) is rather simple: 
we use the backward Euler method of finite differences in order to discretize the parabolic equation 
(fT3l) in time. In each time level we find a new approximation of a solution pair (II, g). First we 
determine a new position of g from the algebraic equation (fl4l) . We remind ourselves that (even in 
the case o is constant) the free boundary function g(r) behaves like rE/q + 0{t 1 ^ 2 ) for r — > + 
(see e.g. or ||28~1 ) and so b(r) = 0(r -1 / 2 ). Hence the convective term b(r)d x Tl becomes 
a dominant part of equation (TT3T) for small values of r. In order to overcome this difficulty we 
use the operator splitting method for successive solving of the convective and diffusion parts of 
equation (fT3l) . Since the diffusion coefficient depends on the solution II itself we make several 
micro-iterates to find a solution of a system of nonlinear algebraic equations. 

Now we present our algorithm in more details. We restrict the spatial domain x E (0, oo) to a 
finite interval of values x E (0, L) where L > is sufficiently large. For practical purposes one 
can take L m 3 as it corresponds to the interval S E (Sf(t)e~ L , Sf(t)) in the original asset price 
variable S. The value Sf(t)e~ L is then a good approximation for the default value S = if L « 3. 
Let us denote by k > the time step, k = T/m, and, by h > the spatial step, h = L/n where 
m,n E N stand for the number of time and space discretization steps, resp. We denote by 11^ an 
approximation of H(xi, Tj), gi rs Q(i~j), V ~ b( T j) wner e Xi = ih, Tj = jk. We approximate the 
value of the volatility a at the node (xj, r/) by finite difference as follows: 



Then for the backward in time Euler finite difference approximation of equation (TT3T) we have 



rr/ ajig'.W) rx ( ( 1 1 / 



Ui)/h^e- x % Tj ) 



w - rp- 



+ (v - l(° j ) 2 ) d x W - \d x ({a j fd x W) + rW = 



k 



(18) 
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and the solution IF = IF(x) is subject to Dirichlet boundary conditions at x = and x = L. We 
set U°(x) = IL(x, 0). Now we decompose the above problem into two parts - a convection part and 
a diffusive part by introducing an auxiliary intermediate step IF - 2 ; 
(Convective part) 

rp'-§ _ rp- 1 . i 

+ ^^n J -5 = o, (19) 

k 

(Diffusive part) 

ni -rp'-i , ,, „ 1 



k - (a j fd x U j - -d x ({aifdJV) + rlF = . (20) 

The idea of the operator splitting technique now consists in comparison the sum of solutions to 
convective and diffusion part to a solution of (TT8l) . Indeed, if d x W ~ d x W~s then it is reasonable 
to assume that IF computed from the system (fT9~l)-(l20"l) is a good approximation of the system (PT8l) . 
The convective part can be approximated by an explicit solution to the transport equation: 

d T fl + b(r)d x fl = for x > 0, r e (tj_i, tj] (21) 

subject to the boundary condition 11(0, r) = —E and initial condition Tl(x,Tj_i) = IF _1 (:r). 
Since the free boundary g(r) = Sf(T — r) must be an increasing function in r and we have 
assumed < q < r we have b{r) = g(r) / g(r) + r — q > and so the in-flowing boundary 
condition 11(0, r) = — E is consistent with the transport equation. Let us denote by B{r) the 
primitive function to b(r), i.e. B{r) = ln^(r) + (r — q)r. Equation (T2T)) can be integrated to 
obtain its explicit solution: 

f l(x r) = { IF"" 1 ^ - B(t) + B{r^ x )) if x - B(r) + Bfa-O > , 
1 ' ; [ —E otherwise. V 

Thus the spatial approximation Ilj 2 can be constructed from the formula 

n i-| _ f n j_1 (6) if^ = Xi -ln^- + ln^_i - (r-g)/c> 0, 
* I — £ otherwise, 

where a linear approximation between discrete values Hp 1 , 2 = 0, 1, n, is being used to com- 
pute the value IF -1 ^ — In Qj + In Qj_ x — (r — q)k). 

The diffusive part can be solved numerically by means of finite differences. Using central finite 
difference for approximation of the derivative d x W we obtain 



TF - TF ^ IF — TF 1 / IF - IF IF - IF 

Hence, the vector of discrete values IF = {Hi, i = 1, 2, n} at the time level j £ {1, 2, m} 
satisfies the tridiagonal system of equations 

afllt! + + 7? nj + ! = Dj"' (24) 
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for i — 1, 2, n, where 



4 s «^n>)-^)> + ±*f 

The initial and boundary conditions at r = and x = 0, L, resp., can be approximated as follows: 



n 



o 



-E for Xi < In (r/g) 
for Xi > In (r/g) 



fori = 0,1,..., n, and IT = -E, IP n = 0. 

Next we proceed by approximation of equation (fl4l) which introduces a nonlinear constraint 
condition between the early exercise boundary function g(r) and the trace of the solution II at the 
boundary x = (S = Sf(t) in the original variable). Taking a finite difference approximation of 
d x Ii at the origin x = we obtain 

= y + ^ 2 ((ni - n$)//», ^ r) HLrJS (26) 



Now, equations (1231) . (|24|) and (|26l) can be written in an abstract form as a system of nonlinear 
equations: 

g> = F(IP, g>) 
IP'^ = T(IP, p>) (27) 
A{U j , ^')IP = U j -^ 

where JF(IP , gP) is the right-hand side of the algebraic equation (l26l) . T(IP, g>) is the transport 
equation solver given by the right-hand side of (|23l) and A = .4. (IF, q*) is a tridiagonal matrix 
with coefficients given by (|25l) . The system (1271) can be approximately solved by means of suc- 
cessive iterates procedure. We define, for j > 1, IF>° = IP -1 , gj>° = gP -1 . Then the (p + l)-th 
approximation of II 3 ' and g> is obtained as a solution to the system: 

n i-i, P +i = T(n J ' p , (28) 
^(n J - p , ^p+ 1 )ip >+1 = ip'-^ +1 . 

Notice that the last equation is a linear tridiagonal equation for the vector IP> P+1 whereas ^' p+1 and 
W~2'P +1 can be directly computed from d26l) and d23l) . resp. Now, if the sequence of approximate 
solutions {(IP' P , £^' p )}21i converges to some limiting value (TP' 00 , g> ,co ) then this limit is a solution 
to a nonlinear system of equations (1271) at the time level j and we can proceed by computing the 
approximate solution the next time level j + 1. 
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Figure 1: a) A comparison of the free boundary function q(t) computed by the iterative algorithm (green 
solid curve) to the integral equation based approximation (dashed red curve); b) free boundary positions 
computed for various mesh sizes; c) a solution profile H(x, r) for t = (blue line), r = T/2 (red curve), 
t = T (green curve); d) contour plot of the function H(x, t). 



4 Numerical approximations of the early exercise boundary 

In this section we focus on numerical experiments based on the iterative scheme described in the 
previous section. The main purpose is to compute the free boundary profile Sf(t) = g(T — t) for 
different (non)linear Black-Scholes models and for various model parameters. We used n = 750 
spatial points and m = 225000 time discretization steps. In average we needed p < 6 micro- 
iterates ((281) in order to solve the nonlinear system (1271) with the precision 10~ 7 . A solution (II, g) 
has been computed by our iterative algorithm for the following basic model parameters: E = 
10, T = 1, r = 0.1, q = 0.05, and a = 0.2. 



4.1 Case of a constant volatility - comparison study 

In our first numerical experiment we make attempt to compare our iterative approximation scheme 
for solving the free boundary problem for an American Call option to known schemes in the case 
when the volatility a > is constant. We compare our solution to the one computed by means of a 
solution to a nonlinear integral equation for g(r) as described in Remark 1 (see also Il2~8ll30"l0 . This 
comparison can be also considered as a benchmark or test example for which we know a solution 
that can be computed by a another justified algorithm. In Fig. [T] part a), we show the function 
g computed by our iterative algorithm for E = 10, T = 1, r = 0.1, q — 0.05, a = 0.2. At the 
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Table 1 : Experimental order of convergence of the iterative algorithm for approximating the free boundary 
position. 



h 


err(L°°) 


eoc(L°°) 


err(L 2 ) 


eoc(L 2 ) 


0.03 


0.5 




0.808 




0.012 


0.215 


0.92 


0.227 


1.39 


0.006 


0.111 


0.96 


0.0836 


1.44 


0.004 


0.0747 


0.97 


0.0462 


1.46 


0.003 


0.0563 


0.98 


0.0303 


1.47 


0.0024 


0.0452 


0.98 


0.0218 


1.48 


0.002 


0.0378 


0.98 


0.0166 


1.48 



expiry T = 1 the value of g(T) was computed as: g(T) = 22.321. The corresponding value g(T) 
computed from the integral equation (see [1281 ) was g(T) = 22.375. The relative error is less than 
0.25%. In the part b) we present 7 approximations of the free boundary function g(r) computed 
for different mesh sizes h (see Tab. Q] for details). The sequence of approximate free boundaries 
gh,h = hi, h 2 , converges monotonically from below to the free boundary function g as h j, 0. 
The next part c) of Fig.[T]depicts various solution profiles of a function Il(x, r). In order to achieve 
a good approximation to equation ((261) we need very accurate approximation of n(x, r) for x close 
to the origin 0. The last part d) of Fig. Q] depicts the contour plot of the function II(a;, r). 

In Tab. CD we present the numerical error analysis for the distance \\gh — g\\ p measured in two 
different norms (L°° and L 2 ) of a computed free boundary position g h corresponding to the mesh 
size h and the solution g computed from the integral equation described in Remark 1 (see also 
[28]). The time step k has been adjusted to the spatial mesh size h in order to satisfy CFL condition 
a 2 k/h 2 ps 1/2. We also computed the experimental order of convergence eoc(L p ) for p = 2, oo. 
Recall that the experimental order of convergence can be defined as the ratio 

eoc(L p ) = ln (H^ ~ gjp) - ln (llgfe t -i ~ q\\p) 
In hi — In hi-i 

It can be interpreted as an exponent a = eoc(L p ) for which we have ||^ — g\\ p = 0(h a ). It 
turns out from Tab. [T]that it is reasonable to make a conjecture that \\g h — g]]^ = O(h) whereas 

\\g h - g\\ 2 = 0(/l 3/2 ) as h -> 0+ 



4.2 Risk Adjusted Pricing Methodology model 

In the next example we computed the position of the free boundary g(r) in the case of the Risk 
Adjusted Pricing Methodology model - a nonlinear Black-Scholes type model derived by Jandacka 
and Sevcovic in [18J. In this model the volatility a is a nonlinear function of the asset price S and 
the second derivative dgV of the option price, and it is given by formula ©. In Fig. [2] we present 
results of numerical approximation of the free boundary position g R (r) = Sf(T — r) in the case 
when the coefficient of transaction costs C = 0.01 is fixed and the risk premium measure R varies 
from R = 5, 15, 40, 70, up to R = 100. We compare the position of the free boundary g R (r) to the 
case when there are no transaction costs and no risk from volatile portfolio, i.e. we compare it with 
the free boundary position q°(t) for the linear Black-Scholes equation (see Fig. [2]). An increase 
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Figure 2: A comparison of the free boundary function g (r) computed for the Risk Adjusted Pricing 
Methodology model. Dashed red curve represents a solution corresponding to R = 0, whereas the green 
curves represent a solution g R (r) for different values of the risk premium coefficients R = 5, 15, 40, 70, 100. 
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Figure 3 : Dependence of the norms 1 1 g R — g° \ \ p (p = oo , 2) of the deviation of the free boundary g = g R {r) 
for the RAPM model on the risk premium coefficient R. 

in the risk premium coefficient R resulted in an increase of the free boundary position as it can be 
expected. 

In Tab. [2] and Fig. |2]we summarize results of comparison of the free boundary position g R for 
various values of the risk premium coefficient to the reference position g = g° computed from the 
Black-Scholes model with constant volatility a = a, i.e. R = 0. The experimental order a p of the 
distance function — g°\\ p = 0(R ap ) has been computed for p = 2, oo, as follows: 

ry — 

p \nRi- In Ri_! 

According to the values presented in Tab.[2]it turns out that the plausible conjecture is that — 
g°\\ p = Oi^R 1 / 3 ) for both norms p = 2 and p = oo. Since the transaction cost coefficient C and 
risk premium measure R enter the expression for the RAPM volatility © only in the product C 2 R 
we can conjecture that \\g R ' c - g '°\\ p = 0(C 2/3 R 1/3 ) as either C -> 0+ or R — > + . 

4.3 Barles and Soner model 

The last example is devoted to the nonlinear Black-Scholes model due to Barles and Soner (see 
A?!). In this model the volatility is given by equation ©. Numerical results are depicted in Fig. [4j 
Choosing a larger value of the risk aversion coefficient a resulted in increase of the free boundary 
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Table 2: Distance || g — g° \\ p (p = 2, oo) of the free boundary position g from the reference free boundary 
position g° and orders and «2 of approximation. 



R 


\\Q H -Q {) \\oo 




\\g H -Q u h 


«2 


1 


0.0601 


- 


0.0241 


- 


2 


0.0754 


0.33 


0.0303 


0.328 


5 


0.102 


0.33 


0.0408 


0.326 


10 


0.128 


0.33 


0.0511 


0.324 


15 


0.145 


0.32 


0.0582 


0.323 


20 


0.16 


0.32 


0.0639 


0.322 


30 


0.182 


0.32 


0.0727 


0.321 


40 


0.2 


0.32 


0.0798 


0.32 


50 


0.214 


0.32 


0.0856 


0.319 


60 


0.227 


0.32 


0.0907 


0.318 


70 


0.239 


0.32 


0.0953 


0.317 


80 


0.249 


0.32 


0.0994 


0.317 


90 


0.259 


0.32 


0.103 


0.316 


100 


0.268 


0.32 


0.107 


0.316 




Figure 4: A comparison of the free boundary function q(t) computed for the Barles and Soner model. 
Dashed red curve represents a solution corresponding to R = 0, whereas the green curves represents a 
solution q(t) for different values of the risk aversion coefficient a = 0.01, 0.07, 0.13, 0.25, 0.35. 
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Table 3: Distance \\g a — g°\\ p (p = 2, oo) of the free boundary position g a from the reference free boundary 
position g° orders and of approximation. 



a 


\\ p a _ pOll 

lie cr 1 1 OO 


Ctrwr, 
■ 


\\Q a ~ £>°|| 2 

lit t: 1 1 Z 


Ct2 


ft 01 


O 1 

u. uo 




ft ftfil 5 




ft r>9 


n 9^ 


O £8 


ft ftQS^ 
U.U70J 


ft £8 
U.Do 


U.UJ 


O 479 


U.U7 


ft 1 84 
U. lo^ 


ft £70 


u.u / 


O f>ft9 


ft 79 


ft 9^9 
U.ZjZ 


ft £Q 


1 
U. 1 


n 7Q^ 


ft 77 
U. / / 


ft 9Q8 


ft 71 9 


n 1 1 


n s^7 

U. OJ / 


ft 89 


ft ^9 


ft 7zL 


0.13 


0.99 


0.86 


0.364 


0.766 


0.15 


1.13 


0.92 


0.409 


0.807 


0.2 


1.52 


1. 


0.529 


0.897 


0.25 


1.97 


1.2 


0.669 


1.05 


0.3 


2.49 


1.3 


0.833 


1.21 


0.35 


3.07 


1.4 


1.03 


1.35 




Figure 5: Dependence of the norms \\g a — g°\\ p (p = oo, 2) of the deviation of the free boundary g = g a {r) 
for the Barles-Soner model on the risk aversion parameter a. 



position g a (r). The position of the early exercise boundary g a (r) has considerably increased in 
comparison to the linear Black-Scholes equation with constant volatility a = a. In contrast to the 
case of constant volatility as well as the RAPM model, there is, at least a numerical evidence (see 
Fig j4] and g a for the largest a = 0.35) that the free boundary profile g a (r) need not be necessarily 
convex. Recall that that convexity of the free boundary profile has been proved analytically by 
Ekstrom et al. and Chen et al. in a recent papers J6l OH dU in the case of a American Put option 
and constant volatility a — a. 

Similarly as in the previous model we also investigated the dependence of the free boundary 
position g = g a (r) on the risk aversion parameter a > 0. In Tab. [3] and Fig. [5] we present results 
of comparison of the free boundary position g a for various values of the risk aversion coefficient a 
to the reference position g = g°. Inspecting values a p of the order of distance || g a — g°\\ p it can be 
conjectured that || g a — g° \\ p = 0(a 2 / 3 ) as a — ► 0. for both norms p = 2 and p = oo. 
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5 Conclusions 



We proposed a new iterative numerical scheme for approximating of the early exercise boundary 
for a class of Black-Scholes equations for pricing American options with a volatility nonlinearly 
depending in the asset prices and the second derivative of the option price. The method con- 
sisted of transformation the free boundary problem for the early exercise boundary position into 
a solution of a nonlinear parabolic equation and a nonlinear algebraic constraint equation. The 
transformed problem has been solved by means of operator splitting iterative technique. We also 
presented results of numerical approximation of the free boundary for several nonlinear Black- 
Scholes equation including, in particular, Barles and Soner model and the Risk adjusted pricing 
methodology model. 

References 

[I] Alobaidi, G. (2004) Laplace transforms and installment options. Math. Models & Methods 
inAppl. Science, 18 (8), 1167-1189. 

[2] Avellaneda, M. & PARAS, A. (1994), Dynamic Hedging Portfolios for Derivative Securi- 
ties in the Presence of Large Transaction Costs. Applied Mathematical Finance, 1, 165-193. 

[3] Barone-Adesi, B. & Whaley, R. E. (1987) Efficient analytic approximations of American 
option values. J. Finance, 42, 301-320. 

[4] Black, F. & Scholes, M. (1973) The pricing of options and corporate liabilities. J. Political 
Economy, 81, 637-654. 

[5] Barles, G. & Soner, H. M. (1998) Option Pricing with transaction costs and a nonlinear 
Black-Scholes equation. Finance Stochast., 2, 369-397. 

[6] Xinfu Chen, Chadam, J., Lishang Jiang, & Weian Zheng (2006) Convexity of the 
Exercise Boundary of the American Put Option on a Zero Dividend Asset, Mathematical Fi- 
nance, to appear. 

[7] Dewynne, J. N., Howison, S. D., Rupf, J. & Wilmott, P. (1993) Some mathematical 
results in the pricing of American options. Euro. J. Appl. Math., 4, 381-398. 

[8] During, B., Fournier, M. & Jungel, A. (2003) High order compact finite difference 
schemes for a nonlinear Black-Scholes equation. Int. J. Appl. Theor. Finance, 7, 767-789. 

[9] EHRHARDT, M. & MlCKENS, R. (2006) Discrete artificial boundary conditions for the Black- 
Scholes equation of American options, submmitted. 

[10] Ekstrom, E. & Tysk, J. (2006) The American put is log-concave in the log-price. Journal 
of Mathematical Analysis and Appl, 314(2), 710-723. 

[II] Ekstrom, E. (2004) Convexity of the optimal stopping boundary for the American put 
option Journal of Mathematical Analysis and Appl., 299(1), 147-156. 



15 



[12] Frey, R. & Patie, P. (2002) Risk Management for Derivatives in Illiquid Markets: A Sim- 
ulation Study. Advances in Finance and Stochastics, Springer, Berlin, 137-159. 

[13] Frey, R. & Stremme, A. (1997) Market Volatility and Feedback Effects from Dynamic 
Hedging. Mathematical Finance, 4, 351-374. 

[14] GESKE, R. & JOHNSON, H. E. (1984) The American put option valued analytically. J. Fi- 
nance, 39, 1511-1524. 

[15] Geske, R. & Roll, R. (1984) On valuing American call options with the Black-Scholes 
European formula. J. Finance, 89, 443-455. 

[16] Hull, J. (1989) Options, Futures and Other Derivative Securities, Prentice Hall. 

[17] HOGGARD, T., Whalley, A.E. & WlLMOTT, P. (1994) Hedging option portfolios in the 
presence of transaction costs. Advances in Futures and Options Research, 7, 21-35. 

[18] JANDACKA, M., & SEVCOVIC, D. (2005) On the risk adjusted pricing methodology based 
valuation of vanilla options and explanation of the volatility smile. Journal of Applied Mathe- 
matics, 3, 235-258. 

[19] Kuske R. A. & Keller, J. B. (1998) Optimal exercise boundary for an American put 
option. Applied Mathematical Finance, 5, 107-1 16. 

[20] Kratka, M. (1998) No Mystery Behind the Smile, Risk, 9, 67-71. 

[21] Kwok, Y. K. (1998) Mathematical Models of Financial Derivatives. Springer- Verlag. 

[22] Leland, H. E. (1985) Option pricing and replication with transaction costs Journal of Fi- 
nance, 40, 1283-1301. 

[23] MacMillan, L. W. (1986) Analytic approximation for the American put option. Adv. in 
Futures Options Res., 1, 119-134. 

[24] Mallier, R. (2002) Evaluating approximations for the American put option. Journal of 
Applied Mathematics, 2, 71-92. 

[25] Mallier, R. & Alobaidi, G. (2004) The American put option close to expiry. Acta Math- 
ematica Univ. Comenianae, 73, 161-174. 

[26] MYNEMI, R. (1992) The pricing of the American option. Annal. Appl. Probab., 2, 1-23. 

[27] Roll, R. (1977) An analytic valuation formula for unprotected American call options on 
stock with known dividends. J. Finan. Economy, 5, 251-258. 

[28] SEVCOVIC, D. (2001) Analysis of the free boundary for the pricing of an American call 
option. Euro. Journal on Applied Mathematics, 12, 25-37. 

[29] SCHONBUCHER, P., & WlLMOTT, P. (2000) The feedback-effect of hedging in illiquid mar- 
kets. SIAM Journal of Applied Mathematics, 61, 232-272. 



16 



[30] Stamicar, R., Sevcovic, D. & Chadam, J. (1999) The early exercise boundary for 
the American put near expiry: numerical approximation. Canad. Appl. Math. Quarterly, 7, 
427-444. 

[31] Jichao Zhao, Corless, R. M., & Davison, M. (2005) Compact finite difference method 
for American option pricing. Journal of Computational and Applied Mathematics, to appear, 
available online 22 August 2006. 



17 



